rey <- function(x, sigma=2){
(x/(sigma^2)) * exp(-x^2/(2*sigma^2))
}
sigma = 2
nsim=5000
X = rep(runif(1),nsim) # inicializar a cadeiaLista 6
MCMC
1 Exercício
Considere a distribuição Rayleigh, com densidade
\[ f(x; \sigma) = \frac{x}{\sigma^2} e^{-\frac{x^2}{2\sigma^2}}, \quad x > 0, \sigma > 0 \]
Use o algoritmo Metropolis-Hastings para gerar da Rayleigh quando \(\sigma=2\). Utilize como candidata \(q(.|x_t)\) o modelo \(\chi^2_{x_t}\). Use um burn-in de 2000 e compute a taxa de aceitação. Compare através de um QQ plot os quantis amostrais e os quantis teóricos \(x_q = F^{-1}(q) = \sigma{-2log(1-q)}^{1/2}, ~0<q<1\).
1.1 Solução
Para ger a Rayleigh usando o algoritmo Metropolis-Hastings, vamos definir a função de densidade Rayleigh e a função de densidade candidata \(\chi^2\) e inicializar a cadeia de algum valor.
Tendo definida a função e condições iniciais, precisamos então simular um número de interações desejadas utilizando como probabilidade de seleção o valor fornecido pelo algoritmo Metropolis-Hastings. A cada iteração, vamos gerar uma amostra da candidata \(\chi^2\) e calcular a razão de aceitação. Se a razão for maior que 1, aceitamos o novo valor; caso contrário, aceitamos com probabilidade igual à razão.
for (i in 2:nsim){
Y = rchisq(1, df=X[i-1]) # amostra da candidata
alpha = min((rey(Y)*dchisq(X[i-1],Y))/(rey(X[i-1])*dchisq(Y,X[i-1])),1)
X[i]= X[i-1] + (Y-X[i-1])*(runif(1)<alpha)
}Vemos pelo gráfico abaixo que a cadeia parece convergir.
plot(3000:5000,X[3000:5000],ty="l",lwd=1,xlab="Iterações",ylab="X")
title(main="Gerar Beta : MCMC - Metropolis-Hastings")Podemos ainda, calcular a taxa de aceitação do algoritmo Metropolis-Hastings. A taxa de aceitação é dada pelo número de vezes que aceitamos uma nova amostra dividido pelo número total de iterações.
acceptance_rate <- round(sum(diff(X) != 0) / (nsim - 1),4)*100Nesse caso a porcentagem de aceitação é de aproximadamente 46.79%, o que é razoável para algoritmos MCMC.
Analisando, como solicitado, o QQ plot dos quantis amostrais e teóricos, vemos que o ajuste é bom, o que indica que a amostra gerada é compatível com a distribuição Rayleigh.
#F^{-1}(q) = \sigma{-2log(1-q)}^{1/2},
qrey <- function(q, sigma=2){
sigma * sqrt(-2 * log(1 - q))}
#Compare através de um QQ plot os quantis amostrais e os quantis teóricos
seq <- seq(0.01, 0.99, by=0.01)
real <- qrey(seq, sigma=2)
amostral <- quantile(X[3000:5000], seq)
qqplot(real, amostral, xlab="Quantis teóricos", ylab="Quantis amostrais")Finalmente, podemos comparar a amostra gerada pelo algoritmo Metropolis-Hastings com a amostra gerada diretamente pela distribuição Rayleigh usando o método de inversão. Para isso, vamos gerar uma amostra de 3000 valores da Rayleigh e comparar os histogramas.
# Função densidade Rayleigh
dRayleigh <- function(x, sigma=2){
ifelse(x > 0, (x/(sigma^2)) * exp(-x^2/(2*sigma^2)), 0)
}
# Geração direta da Rayleigh usando inversão
U <- runif(3000)
X_direct <- sigma * sqrt(-2 * log(1 - U))
# Gráfico comparativo
par(mfrow=c(1,2), mar=c(2,2,1,1))
# Histograma da amostra M-H
hist(X[3000:5000], breaks=50, col="blue", main="Metropolis-Hastings", freq=FALSE, xlim=c(0,10))
curve(dRayleigh(x), col="sienna", lwd=2, add=TRUE)
# Histograma da geração direta
hist(X_direct, breaks=50, col="grey",
main="Geração directa", freq=FALSE, xlim=c(0,10))
curve(dRayleigh(x), col="sienna", lwd=2, add=TRUE)2 Exercício
Considere o modelo normal multivariado \[X_1, X_2, ..., X_p \sim N_p(0,(1-\rho)I + \rho J)\] em que \(I\) é a matriz identidade \(p \times p\) e \(J\) é uma matriz de uns \(p \times p\). Esse é o chamado modelo equicorrelation, uma vez que \(Cor(X_i, X_j) = \rho\) para todo i e j. Note que:
\[ X_i \mid x_{(-i)} \sim \mathcal{N} \left( \frac{(p-1)\rho}{1 + (p-2)\rho} \, \overline{x}_{(-i)}, \; \frac{1 + (p-2)\rho - (p-1)\rho^2}{1 + (p-2)\rho} \right) \]
em que \(x_{(-i)} = (x_1, x_2, ..., x_{i-1}, x_{i+1}, ..., x_n)\) e \(\bar{x}_{(-i)}\) é a média desse vetor. Implemente um amostrador Gibbs usando a condicional acima. Execute o código R para \(p = 5, \rho =0.25\) e verifique graficamente que as marginais são todas N(0,1).
2.1 Solução
Para implementar o amostrador Gibbs, precisamos das condicionais completas de cada variável dada as outras. A condicional foi dada no exercício, e vamos usar essa informação para gerar as amostras.
nsim <- 5000;
X=array(0,dim=c(nsim,5))
rho =.25
p = 5
Std = sqrt((1 + (p-2)*rho - (p-1)*rho^2) / (1 + (p-2)*rho))
X[,]=rnorm(1) # inicializar a cadeia
for(i in 2:nsim){
for (j in 1:p){
Y = X[i-1, -j] # vetor sem o j-ésimo elemento
mean_cond = ( (p-1)*rho / (1 + (p-2)*rho) ) * mean(Y) # média condicional
X[i-1,j] = rnorm(1, mean=mean_cond, sd=Std) # amostra condicional
}
}Abaixo podemos ver o gráfico das marginais, que devem ser N(0,1) para cada variável. Também caulculou-se o p-valor do teste de normalidade de Shapiro-Wilk para cada marginal, que mostrou que todas as marginais são compatíveis com a normalidade.
# Gráfico das marginais com p-valor shapiro
par(mfrow=c(2,3), mar=c(2,2,1,1))
for (j in 1:p) {
# Teste de normalidade Shapiro-Wilk
shapiro_test <- round(shapiro.test(X[,j])$p.value,3)
hist(X[,j], breaks=30, main=paste("X", j,", p = ",shapiro_test), xlab="", ylab="", freq=FALSE)
curve(dnorm(x, mean=0, sd=1), col="red", lwd=2, add=TRUE)
}Por fim, analisando o gráfico de correlação 2 a 2 das variáveis, podemos ver que as variáveis estão correlacionadas conforme esperado pelo modelo (uma correlação positiva, mas fraca). A linha vermelha representa a reta de regressão entre as variáveis.
# crafico de correlação 2 a 2 das pargunais
par(mfrow=c(2,5), mar=c(2,2,1,1))
for (i in 1:(p-1)) {
for (j in (i+1):p) {
plot(X[,i], X[,j], xlab=paste("X", i), ylab=paste("X", j), main=paste("Correlação entre X", i, "e X", j))
abline(lm(X[,j] ~ X[,i]), col="red")
}
}3 Exercicio
Considere o seguinte modelo. \(y_i \sim i.i.d. ~N(\mu, \tau-1)(i = 1, ..., n)\)
em que \(\mu \sim N(0, \beta-1)\) e \(\beta \sim Gamma(2,1)\) e \(\tau \sim Gamma(2,1)\). Implemente um amostrador Gibbs que pode ser usado para estimar \(\mu, \beta, \tau\), em que a amostra observada é dada por:
######################################
## Código R
######################################
y=c(2.34, 2.55, 1.91, 2.42, 2.26, 1.82, 2.04, 1.88, 2.03, 2.28, 2.16, 2.19,
1.89, 1.38, 1.95, 2.50, 1.16, 2.09, 1.46, 0.97, 1.68, 1.89, 2.06, 2.07,
1.72, 1.01, 2.96, 2.84, 2.03, 3.08, 1.80, 2.76, 2.04, 2.05, 1.91, 1.51,
2.00, 0.77, 1.92, 1.54, 2.29, 2.28, 2.20, 1.60, 2.11, 1.36, 1.46, 2.53,
1.89, 1.98)3.1 Solução
Um dos grandes desafios em modelos bayesianos é a amostragem de parâmetros. O amostrador Gibbs é uma técnica poderosa para gerar amostras de distribuições conjuntas complexas, especialmente quando as condicionais completas são conhecidas. Para o modelo dado, as condicionais completas são:
\[ \beta \mid \mu, \tau \sim \text{Gamma} \left( \frac{5}{2},\ 1 + \frac{1}{2} \mu^2 \right) \]
\[ \tau \mid \mu, \beta \sim \text{Gamma} \left( \frac{n}{2} + 2,\ 1 + \frac{1}{2}(n - 1)\, s_y^2 + \frac{1}{2} n \left( \overline{y} - \mu \right)^2 \right) \]
\[
\mu \mid \tau, \beta \sim \mathcal{N} \left( \frac{n \overline{y} \tau}{n \tau + \beta},\ \frac{1}{n \tau + \beta} \right)
\] As funções de amostragem para cada parâmetro são definidas abaixo. A função Gibbs_sampler executa o algoritmo Gibbs, gerando amostras de \(\mu\), \(\beta\) e \(\tau\).
Beta_cond <- function(mu) {
shape1 <- 5/2
shape2 <- 1 + (1/2) * mu^2
rgamma(1, shape=shape1, rate=shape2)
}
Tau_cond <- function(mu, y) {
n <- length(y)
s_y2 <- var(y)
shape1 <- n/2 + 2
shape2 <- 1 + (1/2) * (n - 1) * s_y2 + (1/2) * n * (mean(y) - mu)^2
rgamma(1, shape=shape1, rate=shape2)
}
Mu_cond <- function(tau, beta, y) {
n <- length(y)
mean_y <- mean(y)
mean_cond <- (n * mean_y * tau) / (n * tau + beta)
sd_cond <- sqrt(1 / (n * tau + beta))
rnorm(1, mean=mean_cond, sd=sd_cond)
}
Gibbs_sampler <- function(y, nsim=1000) {
n <- length(y)
mu <- numeric(nsim)
beta <- numeric(nsim)
tau <- numeric(nsim)
# Inicialização
mu[1] <- mean(y)
beta[1] <- 1
tau[1] <- 1
for (i in 2:nsim) {
mu[i] <- Mu_cond(tau[i-1], beta[i-1], y)
beta[i] <- Beta_cond(mu[i])
tau[i] <- Tau_cond(mu[i], y)
}
return(data.frame(mu, beta, tau))
}
parametros <- Gibbs_sampler(y)Os trace-plots abaixo mostram a evolução dos parâmetros \(\mu\), \(\beta\) e \(\tau\) ao longo das iterações do algoritmo Gibbs. Esses gráficos são úteis para diagnosticar a convergência.
# Gráficos para diagnostico
par(mfrow=c(2,2), mar=c(2,2,1,1))
plot(parametros$beta, type="l", main="beta", lwd=1, xlab="Iterações", ylab="Beta")
plot(parametros$tau, type="l", main="tau", lwd=1, xlab="Iterações", ylab="Tau")
plot(parametros$mu, type="l", main="mu", lwd=1, xlab="Iterações", ylab="Mu")Abaixo podemos observar a distribuição de cada um dos parâmetros após a execução do amostrador Gibbs. A partir dessas amostras agora podemos fazer inferências sobre os parâmetros do modelo.
# Histogramas das marginais
par(mfrow=c(2,2), mar=c(2,2,1,1))
hist(parametros$mu, breaks=30, main="Marginal de Mu", xlab="", ylab="", freq=FALSE)
hist(parametros$beta, breaks=30, main="Marginal de Beta", xlab="", ylab="", freq=FALSE)
hist(parametros$tau, breaks=30, main="Marginal de Tau", xlab="", ylab="", freq=FALSE)Podemos observar a trajetória dos parâmetros \(\mu\), \(\beta\) e \(\tau\) em gráficos 2D, onde cada ponto representa uma iteração do algoritmo Gibbs. As linhas conectam os pontos consecutivos, mostrando a evolução dos parâmetros ao longo das iterações.
# Suponha que você já rodou:
# parametros <- Gibbs_sampler(y)
mu <- parametros$mu
beta <- parametros$beta
tau <- parametros$tau
nsim <- length(mu)
par(mfrow=c(1,3), mar=c(4,4,2,1))
##### 1) mu vs beta #####
plot(mu, beta, pch="", xlab=expression(mu), ylab=expression(beta),
main=expression(paste("Trajetória: ", mu, " vs ", beta)))
for (i in 1:(nsim-1)) {
lines(c(mu[i], mu[i]), c(beta[i], beta[i+1]), lwd=1) # movimento vertical
lines(c(mu[i], mu[i+1]), c(beta[i+1], beta[i+1]), lwd=1) # movimento horizontal
}
##### 2) mu vs tau #####
plot(mu, tau, pch="", xlab=expression(mu), ylab=expression(tau),
main=expression(paste("Trajetória: ", mu, " vs ", tau)))
for (i in 1:(nsim-1)) {
lines(c(mu[i], mu[i]), c(tau[i], tau[i+1]), lwd=1)
lines(c(mu[i], mu[i+1]), c(tau[i+1], tau[i+1]), lwd=1)
}
##### 3) beta vs tau #####
plot(beta, tau, pch="", xlab=expression(beta), ylab=expression(tau),
main=expression(paste("Trajetória: ", beta, " vs ", tau)))
for (i in 1:(nsim-1)) {
lines(c(beta[i], beta[i]), c(tau[i], tau[i+1]), lwd=1)
lines(c(beta[i], beta[i+1]), c(tau[i+1], tau[i+1]), lwd=1)
}Por fim, vamos observar as amostras geradas pelo algoritmo Gibbs em um gráfico 3D interativo. Isso nos permite visualizar a relação entre os três parâmetros \(\mu\), \(\beta\) e \(\tau\) de forma completa e entender a região do espaço de parâmetros onde as amostras estão concentradas.
# Remover burn-in (exemplo: 200 primeiros)
burnin <- 200
mu_post <- parametros$mu[(burnin+1):nsim]
beta_post <- parametros$beta[(burnin+1):nsim]
tau_post <- parametros$tau[(burnin+1):nsim]
# install.packages("plotly")
library(plotly)
# Data frame com as amostras
df_post <- data.frame(mu=mu_post, beta=beta_post, tau=tau_post)
# Gráfico interativo 3D
plot_ly(df_post, x = ~mu, y = ~beta, z = ~tau,
type = "scatter3d", mode = "markers",
marker = list(size = 2, color = "blue")) %>%
layout(title = "Amostras Gibbs: mu vs beta vs tau",
scene = list(xaxis = list(title = "mu"),
yaxis = list(title = "beta"),
zaxis = list(title = "tau")))